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ABSTRACT 

A class of complete potential-density basis sets in cylindrical {R, (p, z) 
coordinates is presented. This class is suitable for stability studies of galac- 
tic disks in three dimensions and includes basis sets tailored for disks with 
vertical density profiles that are exponential (e"l^l/^o), Gaussian (e~(*/^o) ) 
or locally isothermal {sech^ {z/Zi^)). The basis sets are non-discrete and 
non-biorthogonal; however, the extra numerical computations required (com- 
pared with discrete biorthogonal sets) are explained and constitute a small 
overhead. The method of construction (and proof of completeness) is sim- 
ple and can be used to construct basis sets for other density distributions 
that are best described in circular or elliptic cylindrical coordinates. When 
combined with a basis set designed for spheroidal systems, the basis sets 
presented here can be used to study the stability of realistic disks embedded 
in massive halos. 



Subject headings: galaxies: kinematics and dynamics — galaxies: structure — instabilities 
— methods: numerical — celestial mechanics: stellar dynamics 
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1. Introduction 

The first step in building a galaxy model is to find or approximate the potential ^p generated 
by a model mass density p. Poisson's equation, 

VV = 47rG'P, (1-1) 

can be solved using basis sets of potential-density pairs {{ipj,Pj) : V^Vj = '^T^Gpj}. A 
given density p is expanded in the basis density functions, 

j 

since Eq. (1.1) is linear, the corresponding potential is 

j 

with the same coefiicients cj. An approximate solution is obtained using finitely many 
terms. 

In this way we can approximate the potentials (and force fields) of mass distributions 
that are not amenable to an exact analytical solution of Eq. (1.1). Moreover, basis expan- 
sions are fundamental to semi- analytical normal mode analyses of stellar systems (Kalnajs 
1977) and to a powerful AT-body simulation technique (Clutton-Brock 1972) that is ideally 
suited to stability studies (Earn & Sellwood 1995). The first step required before imple- 
menting these methods is to find a basis set that is well-suited to the model of interest. 

The eigenfunctions of the Laplacian operator always form a complete biorthogonal 
basis set {e.g., Courant & Hilbert 1953, §6.3; Arfken 1985, §9.4) but they do not always 
converge sufficiently fast to typical potentials relevant for galactic dynamics. It is essential 
that a given model and its normal modes can be represented with a modest number of 
basis functions; otherwise, accurate expansions become prohibitively expensive. 

A variety of useful sets have been derived for flat disks (Clutton-Brock 1972; Kalnajs 
1976; Qian 1992, 1993) and for spheroidal systems (Clutton-Brock 1973; Hernquist & 
Ostriker 1992; Saha 1993; Syer 1995; Robijn & Earn 1996; Zhao 1996). However, the 
literature has apparently been void of sets that are well-suited to disks of finite thickness. 
Such basis sets are needed for studies of realistic three-dimensional (3D) disk galaxy models. 
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There is no evidence that the vertical structure of galactic disks varies significantly with 
radius. The goal of the present paper is therefore to find basis sets suitable for studies of 
3D disks with mass densities of the separable form 



Observations indicate that the distribution of luminous mass in disk galaxies is well- 
approximated by (1.3) with 



results if we demand that the disk be locally isothermal (Spitzer 1942; Binney & Tremaine 
1987, problem 4-25). In (1.4), M is the total mass and a and h are scale lengths. 

New basis sets for 3D disks arc derived in this paper (some of the main results were 
summarized briefly by Earn 1995) and a simple proof of completeness is also provided. 
A derivation of the Laplacian eigenfunctions is reviewed below; the present approach is 
essentially a modification of this procedure, together with a simple trick. Different methods 
for constructing basis sets for 3D disks are discussed by Robijn & Earn (1996). 

2. The Standard Basis 

As a first step, we derive the set of eigenfunctions for the Laplacian operator in cylindrical 
(i?, 0, z) coordinates. Thus we seek solutions of 



p {R, z) = {R) (z) . 



(1.3) 





we obtain from Eq. (2.1) the three ordinary differential equations. 




(2.3a) 



(2.36) 



(2.3c) 



dz'^ 



0, 
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where k and m are separation constants; k is real and can be taken positive while m is an 
integer due to periodicity of the coordinate. These equations (2.3) have the well-known 
solutions 

n (R) = Jm (kR) , (2.4a) 

$ {(P) = e'""^ , (2.46) 
Z{z) =e±('='+^)'^'^ (2.4c) 

where Jm is the cylindrical Bessel function of order m. 

It is easy to see that the eigenvalue A must be negative [multiply Eq. (2.1) by i/j and in- 
tegrate by parts using Green's theorem {e.g., Arfken 1985, eq. 1.98) to get — J (VV')^ dV = 
Xjip"^ dV]. For our purpose of obtaining a basis set, we can impose the further restriction 
that A < — fc^. Thus we may write 2{z) = e*^^, where h is any real number; the eigenvalue 
is then A = — (/c^ -|- /i^). The eigenfunctions are 

iPkmh (i?, 0, z) = Jm (kR) e'"^^ e'^' . (2.5a) 

To each eigenpotential corresponds the density 

Pkmh (i?, 0, z) = Jm {kR) e'^^ e^h' . (2.56) 

By standard theorems {e.g., Courant & Hilbert 1953, §6.3; Arfken 1985, §9.4) these 
functions form a complete, orthogonal basis for the space of square-integrable functions, 
iL^(Il^). Note that to obtain correct units explicitly, a factor GM/L should be appended 
to Jm{kR), where L and M are length and mass scales. 

The standard basis of eigenfunctions (2.5) is very simple, but it is poorly suited for 
numerical work with density distributions that are nearly flat, like disk galaxies. In par- 
ticular, individual basis functions do not satisfy the natural galactic boundary condition, 
— > as i?^ -h 2;^ — > oo. At the very least we need basis functions that decay both as 
i? — > 00 and as 2; — > ±00. 
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3. Basis functions for galactic disks 

We can derive more suitable basis sets by relaxing the conditions we imposed in Eqs. (2.1) 
and (2.2). Rather than demanding eigenf unctions, we shall specify the vertical density 
profile p^{z) in advance and seek solutions of Poisson's equation (1.1) with 

t(;{R,(l),z) =n{R)^{(p)Z{z) , (3.1a) 

Note the difference between this and ordinary separation of variables: here we are free to 
choose the form of p^{z) and we will not have Z{z) oc p^{z), unless we choose p^{z) = e^^^. 

In this case, Eq. (1.1) again separates into three ordinary differential equations, namely 
Eqs. (2.3a) and (2.36) together with 

^-k'Z = P. (z) , (3.2) 

which replaces Eq. (2.3c). 

We are interested in bell-shaped density factors p^{z) that resemble actual vertical 
density profiles of disk galaxies, and we impose the boundary conditions that the potential 
factor Z{z) — > as z ^ ±oo. Many analytical solutions of this form can be found for 



Eq. (3.2). Several simple and useful solutions are listed in P?able 3.1] . The first row gives 
Green's function for Eq. (3.2); when inserted in Eq. (3.1) we obtain the well-known Bessel 
function basis set for fiat disks, used by Toomre (1981) for his fiat disk stability analysis. 
Rows 2-4 give examples that allow us to construct 3D basis sets for realistic disk galaxy 
models. 



The density factors p^{z) in P?able 3.i| do not depend on the radial index /c. As a result. 



the full 3D potential of any disk of the form (1.3) with p^{z) drawn from [Table 3.1] can 



be obtained by a single integration over k, provided the Hankel transform of the radial 
density profile p^(-R) is known analytically: 

/•CX) 

i; {R, z) = p^ (z) / Jo (kR) Z (z) S (k) dk , (3.3) 
Jo 



where S{k) is the Hankel transform of Pj^{R). In this context, the solutions in ]Table 3.1 



have been found previously {e.g., Casertano 1983, Kuijken & Gilmore 1989, Sackett & 
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Table 3.1: Useful solutions of Eq. (3.2) 



Density factor p^{z) Potential factor Z{z) 


5{z) 




£-"1^1 


< 


-fe2^(e-"l^l -fe-'^l^l) , k^a 
^^{l + k\z\)e-''\'\ , k = a 


2 2 

g-a z 


-''"Tkl^ [e'^ erfc( A + az) + e~^^ erfc(^ - az)] 


sech^(az) 


k? + qka 


2 Za 2 2a 

+ e-^"^2Fi(,, f + A; 1 + ^ + J 
2 2a 2 2 a J 



Sparke 1990, Cudderford 1993). The principal contribution of this paper is to show that 
these solutions, and simpler solutions given below, can be used to construct complete 3D 
basis sets. 

If we are willing to let p^{z) itself depend on the radial wave number k then we can 
find somewhat simpler solutions from which to construct complete basis sets ( P?able 372]) . 
Several basis functions for each k are then required to represent a fixed vertical density 
profile such as sech^(a2;). However, this may impose no practical disadvantage for appli- 
cations where the system evolves and/or where only the large-scale normal modes of the 
model are required. 
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Table 3.2: Simpler Solutions with p^{z) depending on k 



Density factor (z) 


Potential factor Z{z) 




-i(l + fc|z|)e-'=l^l 


P e-C^^)' 


-iei/''v^[e'=^ erfc(l/2 + fcz) + 6"'=^ erfc(l/2 - fcz)] 


P sech(fcz) 


kz e^^ — cosh (fcz) log (1 + e^^^) 


P sech^(fcz) 


1 + sinh(fcz) arctan[sinh(fcz)] — ^ cosh(fcz) 



4. New basis sets 

So far we have merely found potential-density pairs; we do not yet have 3D basis sets. 
4.1. Construction 



To form basis sets from the potential-density pairs indicated in P?able 3.1 and Table 3.2 



we make the following observation: the left hand side of Eq. (3.2) is invariant under the 
translation z {z — h) for any /i, so if we replace p^{z) by p^{z — h) and Z{z) by Z{z — h) 
then we have another potential-density pair. This corresponds to shifting the original 
configuration up a distance h. As will be shown in the next subsection, if we let h vary 
from — oo to oo then the functions 



tpkmh (R, (P, z) 
Pkmh {R, (f), z) 



Jm (kR) e'"^^ Z{z-h) 

1 



AtiG 



JmikR) e'^'fp^ (z-h) 



(4.1a) 
(4.16) 



form an iL^-complete basis, where k G IR"^^, m G and Z{z) and p^iz) are taken from 
any row in p?able 3.T1 or [Table 3.2| . Representing a given model amounts to weighting and 
stacking the functions (4.1) along the ^-axis. 
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4.2. Completeness 

It is sufficient to sfiow tfiat any member of tlie standard basis (2.5) can be represented. 
Since our sets differ from tlie standard basis only in tlieir vertical factors p^(^), it is enougli 
to prove tliat for any /i G IR tliere exists a function Wh{h') such tfiat 



Wh {h') [z — h') dh' 



Ahz 



(4.2) 



This Fredholm integral equation of the first kind can be solved explicitly analytically (c/. 
Titchmarsh 1937). Taking the Fourier transform of Eq. (4.2) and using the convolution 
theorem, we have 

1 



(-") h (w) 



^ihz ^iuz dz = 6{u + h) 



(4.3) 



where Wh{u) and p^{u) denote the Fourier transforms of Wfi{h') and p^{h') respectively, 
and 6 is the Dirac delta distribution. It follows immediately that 



Wh {h') 



1 



5(« + /^)^-...'^^ 



1 



Ahh' 



(4.4) 



,27rJ_oo Pz{u) V^pA-h) 
This formula is meaningful provided p^iz) has a Fourier transform that is nowhere zero. 
Table 4.1| gives the Fourier transforms of all the functions in [rable 3.1| and [I'able 3.^ [except 



for arbitrary powers of sech(Q!z); for that transform see Oberhettinger 1973, §1.7.211]. Each 
of these is strictly positive, so our new basis sets are complete in iL^(IR^). 

Table 4.1: Fourier transforms 



Density factor (z) 




5{z) 


1 






e-("^)' 


1 p-u^/ia^ 


sech(az) 


i sech(7ru/2Q!) 


sech^(a2:) 


u csch(7r'u/2a) 
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4.3. Expansion coefficients 

The inner product of two potential-density pairs is defined by minus the interaction po- 
tential energy of the two densities p and p', 

^ ~. j ^^^^ 

Here, the asterisk denotes complex conjugation and the factor L/GM"^ makes the inner 
product dimensionless. All functions are assumed to live in JL^(]El^). In cylindrical coor- 
dinates, 

/oo /"27r /"OO 
/ / p*i;' RdRd(l)dz, (4.6) 
-oo Jo Jo 

where we have set G = M = L = 1 for convenience. 

A basis is a linearly independent set that spans the full vector space of potential-density 
pairs. Thus for any basis set {(V'j, Pj) '■ V^V'j = ^T^Pj} aiid any given density p e E?{W{?) 
there is a set of complex coefficients {cj} such that p = Ylj'^jPj- follows immediately 
that (p I ipf) = Cj (^pj I '(pf), so defining 

Mjf = {pj I V/> , (4.7) 

we have 

where is the jj' element of matrix inverse of M. The potential generated by p is 

■0 = Cjil^j. A basis is said to be biorthogonal if Aijj' is diagonal, i.e., (^pj \ '^y) if and 
only if j = j', and biorthonormal if Mjj' is the identity matrix Sjjr. With a biorthonormal 
basis, the dimensionless expansion coefficients are simply Cj = (p | V'i) = (PilV')- -^W 
finite subset of a basis can be made biorthonormal with the Gram-Schmidt algorithm (e.g., 
Arfken 1985). Biorthogonality is clearly preferable, but without it the only extra effort 
required is the one-time evaluation and inversion of the matrix M-jf (4.7), as emphasized 
by Saha (1993). 

In the present context, where the abstract index j represents the triplet {k, m, h) with 
k and h real, the sum over j refers to a discrete sum and two integrals. 



00 

E 

m=— oo 



-oo ^0 



Ckmh Pkmh dk dh , (4.9a) 



V' = XI / / Ckmh^kmhdkdh. (4.96) 



The matrix M.jj' refers to 



Mkmhk'm'h' = {Pkmh I i'k'm'h') , (4-10) 

and biorthonormality means strictly that 

Mkmhk'm'h' = S{k- k') 5^^, 6{h- h') . (4.11) 
Using the Bessel closure relation, 

/•OO 1 

/ Jm (kR) Jm (k'R) RdR = - d{k-k') , (4.12) 
Jo k 

{e.g., Arfken 1985, Eq. 11.59) and the exponential identity, 

/ e-^'^'V™'<^d(^ = 27rW, (4.13) 
we find for the present class of basis sets (4.1), 

Mkmhk'm'h' = S{k- k') 5mm' Mkhh' , (4-14) 

where 

Mkhh' = 1^ pI {z-h)Z{z- h') dz . (4.15) 
The expansion coefficients are therefore 

°o roo poo 

Ckmh = J2 Mkmhk'm'h' iPl'I'k'm'h') dk'dh' 

/OO 

Mkih'ipli^kmh') dh', (4.16) 
-OO 

where Al""*^ is the matrix inverse of M.. To find we must solve the integral equation 

/OO 

MkL"^kh"h' dh" ^5{h-h') . (4.17) 
-OO 

We discuss how to solve this equation numerically in the next section. 
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5. Discretization 

We have given several crucial formulae that involve integrals over the radial and vertical 
basis indeces, k and h [(4.9), (4.16), (4.17)]. Since our goal is to use as few basis functions 
as possible to obtain a given accuracy, these integrations must be reduced to computations 
with a small, finite number of k and h values. 

5.1. Numerical quadrature 

The continuous indeces k and h can be replaced with discrete indeces n and £ via prescrip- 
tions of the form 

N 

f{k)dk ^Y^w^Jikn), (5.1a) 

/oo L 
9{h) dh -^Y.w^^g{h^) , (5.16) 

£=0 

where the sums approximate the integrals to some given order in 1/N and 1/L. The 
simplest approach is to choose a set of evenly spaced abscissae, kn = nAk and hi = 
{£ — L/2)Ah, and to use a classical formula such as Simpson's rule to define the weights 
{e.g., Abramowitz & Stegun 1972, §25.4.6; Press et al. 1992, Eq. 4.1.13). In the case of 
the axisymmetric (m = 0) integrals over k in Eq. (4.9), we can do better by employing a 
Gaussian quadrature rule {e.g., Press et al. 1992, §4.5) because we know the initial radial 
profile. 

5.2. Solving the integral equation for 

In any application, our first numerical task is to solve Eq. (4.17) for To reduce 

notational complexity, let us define 

^ll' ^ Mk^heh, , Bl\> ^ M]^Xh, ■ (5-2) 
We wish to solve for the matrices B^, n = 0, . . . , N. Using (5.1), Eq. (4.17) becomes 

L 

J2^l"BiMl"l' =0, iy^i', (5.3a) 
£"=0 
L 

w\Y.w\,Sl.J^l.i =1, (5.36) 
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where the second hne is obtained by integrating Eq. (4.17) in a small interval about h£ for 
£ = £'. Putting these two pieces together we have 

^ wt„B^^„A1„i, = —Sw . (5.4) 
£"=0 '^i 

Thus is the ordinary matrix inverse of A^W, where W = dia,g{wQ , . . . , w^l) . Each 
matrix needs to be computed only once for a given basis set, so the work involved is 
always negligible. 

5.3. Discrete expansion formulae 

We can now rewrite our formulae for the expansion coeflBcients and expanded potential in 
a completely discrete manner. Eq. (4.16) becomes 

Cnmi = ^ W^i B^^, {p I Ipnml') ■ (5-5) 

e' 

Eq. (4.96) becomes 

mi n 

and a similar formula replaces Eq. (4.9a). Note that in the case of a distribution of point 
particles, p = Mi 5{x — Xi), the inner product in (5.5) also becomes a discrete sum, 

(P I ^nmt) = -^Mi Ipnrne (Xi) ■ (5.7) 



6. Discussion 

A similar approach can be used to derive basis sets for 3D elliptic disks. In fact, all we 
need to do is replace the factor Jm{kr) e^'^'^ of our basis sets with the eigenfunctions of 
the Laplacian in elliptic coordinates (Mathieu functions, e.g., Abramowitz & Stegun 1972, 
chapter 20). Such basis sets should be useful for stability studies of realistically thickenned 
versions of fiat elliptic disks {e.g., Evans & de Zeeuw 1992). 

For stability work, it is essential to be able to represent the normal modes of the 
system with a small number of terms of the chosen basis expansion. There is no guarantee 
that a basis selected for its ability to represent the unperturbed model will be ideal for 
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representing its normal modes. It is prudent, therefore, to repeat stability calculations 
with several different basis sets. Robijn & Earn (1996) discuss alternatives to the present 
class. 

Both in stability analyses and A^-body experiments, only the basis potentials are used 
for the principal computations: in Eqs. (5.5), (5.6) and (5.7) the basis densities do not 
appear. For this reason, it is unfortunate that in Tables 3.1 and 3.2 the potential functions 
are more complicated than the associated densities. Nevertheless, they are easy to im- 
plement because the special functions they involve are available in all standard numerical 
libraries. Moreover, all the potentials and densities discussed in this paper are separable 
so they can be interpolated efficiently and accurately from one-dimensional tables. The 
basis sets provided here should, therefore, be useful for stability studies and modeling of 
disk galaxies. 

I am grateful to Tim de Zeeuw, Ed Doolittle, and Alar Toomre for helpful discussions 
and correspondence, and to Frank Robijn both for discussions and detailed comments on 
a preliminary version of the manuscript. I also thank the referee, Agris Kalnajs, for useful 
comments. I was supported by a Lady Davis Postdoctoral Fellowship. 
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